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, The competition between long-range and short-range interactions among holes moving in an 

■ antiferromagnet (AF), is studied within a model derived from the spin density wave picture of 
^S) ' layered transition metal oxides. A novel numerical approach is developed which allows one to solve 

the problem at finite hole densities in very large systems (of order hundreds of lattice spacings), 

■ albeit in a quasiclassical limit, and to correctly incorporate the long-range part of the Coulomb 
interaction. The focus is on the problem of charge ordering and the charge phase diagram: at low 
temperatures four different phases are found, depending on the strength of the magnetic (dipolar) 
interaction generated by the spin-wave exchange, and the density of holes. The four phases are 
the Wigner crystal, diagonal stripes, a grid phase (horizontal-vertical stripe loops) and a glassy- 
clumped phase. In the presence of both in-plane and out-of-plane charged impurities the stripe 

' ordering is suppressed, although finite stripe segments persist. At finite temperatures multiscale 

, (intermittency) dynamics is found, reminiscent of that in glasses. The dynamics of stripe melting 

QQ ■ and its implications for experiments is discussed. 

cn 

d • I. INTRODUCTION 

, Charge ordering in layered transition metal oxides has recently attracted a significant research interest, due to. its 

^ 1 possible relation to the mechanism of high temperature superconductivity in doped cupratesEJ and bismuthates.0 In 
Q ' particular, stripe-like ordering, which involves holes ordered into linear arrays, separated by an antiferromagnetically 
Q : (AF) ordered electronic backgrouncL have been discussed as a candidate for the explanation of pseudogap effects 
^ in underdoped cuprate compoundsJJ In addition, the formation of domain walls has been discussed in terms of 
• i-H , the proximity to phase separation. a Quite generally, phase separation on mesoscopic and even macroscopic scales 
' is potentially relevant forj-any strongly correlated orgaiH|C and inorganic electronic system^-nincluding systems with 
H ^ spin-density-wave (SDW)j3 charge-density-wave (CDW)P and Jahn- Teller broken-symmetryla ground states. 

On the experimental side, mesoscopic (nanoscale) phase separation has been observed in many compounds. In the 
case of La2-a;Sra;Ni04+y stripes have been obser\£ed both using nuclear magnetic resonance (NMR) methods and more 
directly, using high-resolution electron diffractionlZl. In addition, stripes have also been identified in Lai-ajCa^^MnOs for 
specific commensurate values of dopingp. In cuprates static stripe order has been observed in Lavr&-2:Sra;Ndo.4Cu04 
in both elastic and inelastic neutron scattering experimentsEI, and x-ray diffraction experimentsEj. There are also 
evidences that stripes exist in some form in high- Tc compounds. In the oxygen doped La2Cu04+Jlj stripeaJaave been 
observed using the nuclear magnetic resonance (NMR) techniques. Magnetic susceptibility measurementsllj, nuclear 
quadrupole resonancetj (NQR) and muon spin resonancetil all indicate formation of domains in La2-xSr3;Cu04 and 
recent inelastic neutron scattering (INSi experiments in La2-a;Sra;Cu04 and YBa2Cu307_5 superconductors yield 
results consistent with stripe format ionEjllJ, although the width of the INS lines in, e.g., YBCO materials is large, 
which may suggest dynamic charge ordering. 

On the theoretical side, stripes have been proposed by several research groups. Since in strongly correlated systems, 
such as cuprate superconductors, electrons exhibit a strong on-site repulsion, numerous studes have been devoted to 
the Hubbard and t-J models. ItJias been shown that a mean-field treatment of the Hubbard model yields a stripe 
phase as a locally stable solutionis. Many other studies view the stripes as an outcome of the competition between 
kinetic energy of holes and exchange energy of spins alone and frequently neglect the role of the long-range part 



1 



of the Coulomb interactionJia'EJ and only recently, an|-a|ttempt to incorporate the long-range forces into the mean- 
field approach to the Hubbard model has been madeO Another point of view emphasizes the intria^jie instability 
of a strongly correlated electronic system towards a phase separation as a necessary starting point £3'E3 Then it is 
assumed that such an instability is prevented by the long-range Coulomb forces. Therefore, the competition between 
this instability, whose existence in the physical range of parameters of the realistic models is yet to be proven, and 
Coulomb repulsion gives rise to a stripe phase. Thus, these two approaches agree on the importance of the correlations 
but disagree on the role of long-range forces. More recently, it has been shown that phase separation is indeed a very 
common phenomenon close to quantum critical points.Eil 

One would expect that the existence of stripes in the widely studied "minimal" t — J or Hubbard models can 
be either proven or disproven by some unbiased numerical technique. Unfortunately, numerically the stability of the 
stripe phase has been established less clearly. Numerical studies of the t-J model are presenting conflicting conclusions 
as to the existance of stripe phases ii>-tJa§ ground state of this "basic" strongly correlated model which might be the 
result of the strong finite-size effectscJcj. For example, even a Monte Carlo simulation of the doped Ising model, 
without the long-range forces, yields holes ordered into loops, rather than into geometric arrays.cZl In factj-with an 
exception of the recent Density Matrix Renormalization Group (DMRG) simulations of White and Scalapinoc3, which 
have found stripe formation in relatively large t — J clusters, no microscopic calculation to date has shown that the 
stripes are a stable entity. Most importantly, no stripe formation in a system with long-range interaction has been 
studied in a direct simulation. In addition, the sizes of the clusters available with modern day computers for solving 
quantum models of spins and holes are still too small to study role of the the long-range Coulomb interaction and be 
free from significant finite-size effects. 

In this situation we propose a different strategy: one can study a quasiclassical limit of the quantum problem 
of holes in an AF spin environment analytically and incorporate all essential correlations in an effective hole-hole 
interaction. In this case the AF background is effectively integrated out, and the focus is on the charge subsystem. 
Then the motion of "classical" holes at finite density, interacting via an effective magnetic interaction and in the 
presence of long-range Coulomb forces, can be studied numerically in much larger systems. In other words, in this 
paper we combine analytical and numerical approaches to study the charge ordering in transition metal oxides. 

Our numerical-approach is based on the spin density wave (SDW) picture of Schrieffer_Wen and Zhang,E3 which 
is closely relatedtJ to the semiclassical approach to the t-J model by Shraiman and SiggiatHl in which the interaction 
between doped holes stems from the spiral distortion of the local Neel vector near a hole. As shown below, in 
the quasiclassical limit, the problem can be solved using classical Monte Carlo (MC) or molecular dynamics (MD) 
methods. In a systematic numerical study we explore the interplay between long-range Coulomb interaction and 
short (or intermediate)-range AF interactions of dipolar nature, which we take to have both isotropic and anisotropic 
components (depending on the lattice structure). 

Our main results can be summarized as follows: in the absence of disorder we find four phases depending on the 
density of holes and the characteristic AF energy scales: (z) a Wigner crystal, {ii} diagonal (glassy) stripes, (m) 
a geometric phase, characterized by horizontal- vertical stripes or checkerboard (grid), and {iv) a "clumped" phase 
(phase separation). In our study the stripe- like phases emerge as a kind of melting of the Wigner crystal phase, hence 
the long-range Coulomb interaction is a necessary ingredient for their occurrence. In the geometric phase the stripes, 
resulting from the competition of the short-range and long-range interactions, are characterized by a particular AF 
dipolar alignment. The patterns are very stable, showing large "string tension," while the motion of holes within a 
stripe is much softer. If one takes into account the kinetic energy of the holes along the stiioes one is lead to the 
concept of a quantum liquid crystal as proposed recently by Emery, Fradkin and KivelsonH On the other hand, 
the ground state of the geometric phase is not well defined in that there are many geometric phases with very low 
energies, comparable to that of the ground state, implying a rugged energy landscape. We find that, on lowering the 
temperature, the geometric hole ordering is characterized by occurence of secondary defects in the structure. At higher 
temperatures we find that the dynamic hole ordering is characterized by temporally intermittent pattern formation 
(i.e., spatio-temporal intermittency). Finally, we find that a sufficient concentration of randomly placed impurities 
destroys the geometric hole pattern, although, regardless of the impurity type, stripe segments are preserved. 

The paper is organized as follows: in the next section we present a review of the theoretical model and the 
computational methods we use. In section IH we present our numerical results and in particular we present the phase 
diagram showing how the obtained phases emerge as a function of doping and interaction strengths. Finally, in section 
IV we summarize our conclusions and experimental implications. 
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II. MODEL 



We begin with the spin density wave (SDW) picture of layered transition metal oxides. This pictupi.iias been 
very successful in describing the stoichiometric insulating AF phase of these systems at low temperaturesoo In this 
picture the electrons move with hopping energy t in the self-consistent staggered field of its spin, as described by, e.g., 
the Hubbard model (the calculation is presented in detail in Appendix ^ : 

H = -tJ2i4cj + h.c.) + C/ ^ n,^]n,^^. (1) 

Because, the translational symmetry of the system is broken, the electronic band is split into upper and lower Hubbard 
bands.Ej On performing a Bogoliubov transformation, one defines the valence (hk^a) and conduction band {pk,a) 
operators, respectively: 

Pk,a = UkCk,a + OlVkCk+Q.a (2) 
hk,a = UkCk,a - aVkCk+Q.a, (3) 

where Ufc, Vk are the Bogolioubov weights and a is the sublattice index. The upper and lower Hubbard bands are 
separated by the Mott-Hubbard gap, A = US/2, where S is the expectation value of the staggered field Sz: 

{SM)) = -'^'^Uk+q-QVk{hk+q-Q^chl,^), (4) 
k.a 

calculated at momentum transfer q = Q. At half filling the lower band is filled and the upper band iSpempty. This 
picture is consistent with the angle resolved photoemission data in the layered AF insulator Sr2Cu02Cl2.tJ On doping 
the system with holes with planar density ct^, at low temjpGjratures, T <C A./kB, the low frequency physics reflects 
purely the lower Hubbard band (LHB). It has been shownEJ that, regardless of the band structure, the LHB has a 
maximum at four wavevectors = (±1, ±l)7r/ (2a), where a is the lattice spacing, and therefore the long wavelength 
theory of the problem can be studied by assuming the momentum of the holes to be close to these points. 

Then the two hole interaction Hamiltonian can be separated into the longitudinal and transverse parts (Hz and 
Hxy, respectively), whose Fourier transform, for quasiparticle momenta near ki, is equal to:ca 

(r)= [Azala'2 - A,y {a+a^ + a^(J+)] S{r) - 



di ■ da _ ^ (di ■ r) (da ■ r) 



{cj+a^ + (y+) , (5) 



where r — |ri — ra] is a relative hole-hole distance, is a coordinate of a hole in units of a, crf^^'' = cJ,(t^^^^''c/3 is a 

spin-density operator, with cr^^^-' Pauly matrices, d^ is a unity vector in the direction of the dipole moment of the 
hole. In the SDW formalism the interaction strengths A^, A^y and B^y ^ Axy/^n are all of order Hubbard U . The 
interaction (^) is clearly rotationally invariant and valid for r > a, while for r — s- it yields an unphysical divergence 
of the (attractive) dipolar interaction. We observe that this form of the Hamiltonian is not particular to the SDW 
theory of the Hubbard model, but stresses the fact that a mobile carrier in an antiferromagnet produces a dipolar 
distortion of the magnetic background. We demonstrate this explicitly in Appendix ^ where we show that the t-J 
model has exactly the same type of interaction terms. In other words, at finite density the holes interact via two 
different mechanisms: a uniform short-range attractive force due to AF bond-breaking and a long-range magnetic 
dipolar interaction (contained by Eq. (H), see Ref. ^). The latter term is due to the long range spiral distortion 
of the Affj-background, which is a consequence of quasiparticles interacting with soft (Goldstonc) modes of the spin 
system.cj'cJ The magnetic dipole moment associated with each hole is due to the coherent hopping of holes between 
different sublattices and scales with the AF magnetic energy. This implies that the quantum effects associated with 
hole kinetic energy can be neglected, which is correct in the limit t <C J, believed to be valid in nickelates. This is 
also why the hole-hole interaction obtained ia, the weak coupling SDW picture is equivalents to that in the effective 
Hamiltonian found by Shraiman and SiggiatO, based on the t — J model, where the dipolar interaction is obtained 
using semiclassical analysis of the spin part of the model, as well as symmetry considerations. It is also possible 
to prove, using Ward identities, that the remaining spin part of the problerp-js equivalent to the two dimensional 
(2D) non-linear a model in the long wavelength limilCJ. It has been arguecO that at physical values t/ J ^ 1 all 
coupling strengths {A^, A^y and B^y), and therefore the hole- hole interaction, will be renormalized to the value of 
super-exchange constant J. 
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In pure two dimensions at finite T the system is magnetically disordered, characterized by a finite magnetic cor- 
relation length, ^ (see Ref. and the range of the dipolar interaction between the holes, mediated by the AF 
background, is also of order ^. In fact, even at T = and finite hole concentration the correlation lenght is restricted 
and the dipolar intejcaction is effectively short ranged. 

It has been notedE3 that the dipolar twist of the magnetic background would imply local time-reversal symmetry 
breaking which puts some restrictions on the applicability of the picture to the real systems where such symmetry 
breaking has not been experimentally observed. Indeed, the time-reversal symmetry is already broken in the Neel 
state, as well as in any magnetically ordered state. However, in two-dimensional spin systems at finite temperature 
all symmetries are restored and we expect this to be true also for the hole-doped case studied here. 

Besides the AF interactions the holes also experience the long-range Coulomb interaction. This is clear if we 
consider that Vg = rg / ag (where tq is the mean inter-particle distance and Og is the Bohr radius) is very large in the 
underdoped systems (r^ « 8). In other words, in systems with a small density of holes the screening, which is due to 
the density fluctuations, is very weak and one must take the Coulomb interaction into account. 

Finally, each hole carries a spin degree of freedom as well. However inspection of Eq. (|^) reveals that the overall 
spin energy is minimized in the spin anti-symmetric channel. Hence we neglect the spin-symmetric channel and thus 
in our approach, we only consider the charge channel with an effective (magnetic in origin) interaction between two 
holes, I and 2, of the form (see Eq. (||)) 

„2 

V{y) = ^ - Ae-'''" - B cos{2e ~ (j}i - (j)2)e-'''^ . (6) 
er 

where we have asssumed that r can be relaxed from a crystal lattice position to an arbitrary (continuous) value. We 
return to this point later. 

In Eq. (|^) q is the hole charge, e is the dielectric constant (which we assume to be of order I), 61 is the angle made 
between r and a fixed axis and (j)i^2 are the angles of the magnetic dipoles relative to the same fixed axis. A is the 
strength of the uniform (short-range) interaction and B is the strength of the magnetic dipolar interaction, which we 
will assume to be independent adjustable parameters, which, in real materials, should be of order ~ I eV. Note that we 
have introduced B ~ B^y/l^, where I is some appropriate average length, a < ^ < ^, in order to avoid the unphysical 
divergence of the dipolar part of the interaction in Eq. (|^), while keeping the necessary symmetry of the interaction. 
Moreover, in the SDW picture, at low doping, and (j)2 are restricted to the angles {2n + I)7r/4, the four angles 
determined by the vectors k^. In addition, the dipole moment vectors are also restricted to ki, which justifies our use 
of Eq. (^ where we have assumed a fixed size of the dipole moment for each hole. However, at larger doping levels 
these angles can be relaxed to arbitrary values, provided the interaction is short-ranged. Of course, we have verified 
by an explicit calculation that restricting dipole angles to the discrete values does not qualitatively change our results, 
presented in the next section. It is interesting, to note that the hole-hole interaction in the form almost identical to 
Eq. (^ has been obtained by Aharony et alx3 for the static holes residing on Cu-0 bonds within the framework of 
a classical model of an AF diluted by ferromagnetic bonds. In this case the value of the coupling constant B is also 
restricted by a few J. Indeed, starting from the insulating phase of cuprates, the holes are injected into the Cu02 
planes at high temperatures during the sample preparation. The hole distribution in this case is annealed (instead 
of quenched as proposed in Ref. |39|) since the holes have enough phase space for interactions. As the temperature is 
lowered the holes can adjust themselves to V{r) and form the structures we discuss below. 

Quite generally, one can think of an AF as an active media generating long range dipolar forces in response on 
some local distortion. Therefore, the interaction ^(r), Eq. (^) is of more general significance than just a result of the 
SDW picture, and the study of the system of particles interacting via V{r) is of wider interest. 

In general, the many-body problem of holes in an AF background is extremely complicated, involving many- 
particle interaction terms. However, at low densities, where the average distance between holes is comparable to the 
AF correlation length, it is reasonable to assume that the interaction of any two holes is weakly perturbed by other 
holes, and the total potential energy can be expressed in terms of two-particle energies, provided the Af correlation 
length is replaced by an effective correlation length, which, to avoid clutter, we also denote as ^. We therefore study 
the system of "classical" particles in a computational box of size Lx x Ly, interacting via a potential 

<i,j> 

where V{r) is given by (^. However, we emphasize that our approach is not a self-consistent one in the sense that 
the true interaction must include many particle terms (omitted here), which stem from the fact that the SDW state 
is altered due to the charge ordering. The self-consistent approach to charge ordering will be presented elsewhere. In 
addition, superconducting fluctuations have been neglected. Moreover, the kinetic energy of the holes may lead to a 
quantum melting of the phases discussed here. 
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FIG. 1. Top: Electrostatic interaction in a system with periodic boundary conditions: an effective interaction between any 
two charged particles in the computational box (central rectangle), such as the two marked by a thick vector line, involves the 
interaction with all of the particle images, marked by thin wave vector lines. Bottom: The Lekner potential (which accounts 
for the periodic boundary conditions) in comparison with the bare Coulomb potential, within a computational box of size 
L. As expected, at r — » the two potentials coincide (up to an additive constant) and hence their difference can be readily 
approximated by a low order polynomial. 

The long-range character of the Coulomb interaction requires further consideration: in order to perform calculations 
at finite density, as required by the dipolar interaction we consider, we usually perform calculations with periodic 
boundary conditions. The ability to handle long-raage Coulomb interactions in rectangular periodic media has been 
enhanced recently in the area of molecular physics. CJ Assuming a computational cell of arbitrary geometry and cyclic 
boundary conditions it is possible to sum interactions .of. particles with all of their images residing in cells obtained 
by translation from the original computational cell,Eilc3 yielding an effective interaction which is periodic in the 
computational cell used (see Fig. |l]). On making integral transformations. Coulomb interactions are computed by 
summing over fast-convergent Bessel functions with great accuracy (see Ref. ^ for a detailed study of the Lekner 
summation technique). The computational efficiency is further enhanced by tabulating the effective interaction. This 
is possible since the difference between the obtained effective interaction and the Coulomb interaction is a well behaved 
function which can be easily calculated using polynomial interpolation. 

At the beginning of each simulation we place the holes at random and assign to each hole a magnetic dipole moment 
of constant size, but random direction. We study the energy landscape and the dynamics of the system aosing three 
different methods: Monte Carlo (MC), Langevin molecular dynamics (MD) and a hybrid MC-MD method-c3 All three 
methods yield essentially the same results. Both MC and MD methods are well known in the literaturecj and hence 
we only review details pertinent to our calculation. In the MC method we use the standard Metropolis algorithm for 
the acceptance of hole configurations. For the Langevin MD method the dynamics of the system is determined by the 
forces, obtained from Eq. (|j), with a noise term, rji, for each particle which satisfies the usual fluctuation-dissipation 
condition (77^77^) = 27TjTSij, where i, j correspond to all possible degrees of freedom and 7 is a damping term.c3 Since 
the system exhibits several phases (see Fig. p|) for some values of the input parameters, its ground state is not always 
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well defined and a numerically obtained low temperature state may, in fact, depend on the initial and boundary 
conditions. Hence, in order to rapidly reach a hole configuration with the lowest global minimum energy we perform 
simulated annealing from high temperatures. 

The hybrid MC method includes elements of both the MC and MD methods: the hole configuration is again deter- 
mined using the standard Metropolis algorithm, but here a new configuration is obtained by letting the system evolve 
through a classical MD calculation over a certain time period. Note that, in principle, in classical MD calculations the 
energy is a conserved quantity, hence every step should in principle be accepted. In reality the MD method introduces 
errors which typically slightly lower the system enepgy, just as required by the Metropolis algorithm. Hence this 
method yields extremely high MC acceptance ratios.EJ 

III. RESULTS 

In this section we discuss the results of our numerical calculations. We first present the obtained ordered phases 
and the phase diagram of the system and discuss its implications. We then show the hole ordering in the presence 
of disorder. Finally we show the properties of the system as a function of T and in particular the dynamics of the 
observed (stripe) pattern formation. 

Before we begin with the presentation of our results we address the input parameters of the model, namely the 
hole density, ct^ (or the doping level n) , the strength of the isotropic and dipolar part of the AF interaction, A and B 
respectively, the AF correlation length, ^, the temperature T and the concentration of impurities, Ci, for the systems 
with static point disorder. We define the doping level n as the hole density measured in units of cuprate lattice 
spacing; thus n — 1% corresponds to 1 hole per 100 a^, where a « 3.8 A. 

The input parameters are not necessarily independent of each other, e.g^ A and B should be proportional to each 
other, with A ^ U when t ^ U and A « it'^/U when U ^ t (see Ref. However, since the range of the bond 
breaking and dipolar interactions is vastly different, it is reasonable to treat A, B and ^ as independent parameters. 
Clearly both A and B should be of order of the Hubbard U in the SDW approach and of the order of J in the strong 
coupling limit, and the correlation length is of order 1 to 10 lattice spacings in cuprates. 

We begin with the low temperature properties (ground state) of the system as a function of B. The relevant order 
parameter for charge ordering is the Fourier transform of the hole density: 

1 ^ 

P^'i^ = nT.^'^""' ^ (8) 

1=1 

where is the position of the i*'* hole and N is the total number of holes. A peak in p(q) at some wave- vector q = K 
indicates ordering. Returning to the SDW picture, presented in Sec. ^ we recall that the hole density is given by 

From the definition of the staggered magnetization, Eq. (^), it is immediately clear that {Sz{q)) oc Sg^K+Q, i-e-, a peak 
in pk at K leads to a magnetic peak at Q + K and by symmetry at Q ~ 

Since the interaction due to the second term in Eq. (^ (isotropic attractive interaction) is extremely short-range 
(in fact in an infinite system it is a (5 function), it is initially reasonable to set A — and explore the behavior of 
the system as a function of B. We return later to the role of A. As explained in Sec. || low T properties have been 
obtained by annealing the system from some high temperature (T ~ 5000K) down to temperatures of order IK. la 
the extreme case B — we find a Wigner crystal with small distortions, to be the state of lowest energy, as expectedo 
(see Fig. |a). 
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FIG. 2. Low T states of the hole vacancy system, (a) For 5 = the holes (circles) form the Wigner crystal and (b) for 
_B — » oo (an unphysical case) they form a "clump" pattern, a characteristic of the "mesoscopic phase separation." The dipole 
orientation (shown by the segments, originating from the circle centers) indicates finite magnetization at each star cluster. In 
both panels the doping level is n = 15%. 



The small distortion of the Wigner crystal structure is due to the periodicity, which introduces a small spatial 
anisotropy into the system due to the shape of the computational box. Indeed, upon setting, e.g., Ly/L^ = •\/3/2, 
we obtain a perfect Wigner crystal to be the ground state. Another extreme case is when B oo. In this case the 
AF dipolar interaction dominates over the average Coulomb interaction; one then finds star shaped clumps of holes, 
similar in shape to those found in Ref. 45, which can, at sufficiently high density, form a geometric structure (e.g., a 
Wigner crystal of clumps). We note that this case is rather unphysical, as macroscopic phase separation is inconsistent 
with our initial assumption of the two-body dipolar interaction being independent of the many-hole effects. On the 
other hand, ordering of niaiiroscopically hole rich regions is in agreement with the conjecture that all ground states 
are geometrically ordered O 

On increasing _B > 0, at fixed density, the Wigner crystal becomes unstable and a new phase with diagonal stripes 



is formed, as shown in Fig. 2a of Ref. 47. The main characteristic of this phase is a ferro-magnetic ordering of the AF 
dipoles. The situation here is very similar to that observed in La2-xSr:rNi04+y (see Ref. ^). Note that such a state 
with a dipole ordering appeares to violate time reversal symmetry. On the other hand the true ground state in this 
case also involves hole ordering, with holes aligned in stripes either perpendicular or parallel to the dipole orientation. 
However, the interstripe distance, in this case, is close to that between holes within a stripe and hence a simulation 
inevitably yields a "glassy" state, with many defects. This is reflected in the shape of pk, which shows broad peaks 
(see Fig. 2b in Ref. E^, indicating an average interstripe distance. 
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FIG. 3. Low T state for larger finite values of B, B = 4el/, at n — 15% doping level. The holes (circles) form a gridl 
(panel (a)) with dipoles orientated along line segments and a dipole rotation at grid intersections (panpl-|(b)). Panel (c) shows 
a contour plot of the average charge density pk in arbitrary units, indicating "perfect" geometric order JIZI Even after averaging 
over many solutions in this case the charge peaks are much sharper than those found in the ferro-dipolar phase (see Fig. 2b in 
Ref. 1^). 

As shown in Fig. |^a, at larger values of B a linear stripe is formed, which, with increasing density tends to close 
into a loop. More importantly, the loop formation is accompanied by magnetic dipole orientation along the straight 
portion of a loop with gradual rotation by tt/2 at each corner-ea Due to the rotation of dipoles at corners the loops 
interact, and eventually form the checkerboard (grid) pattern.EJ The size of the distance between holes within a line 
is determined by the ratio of B (or the sum of A and B, for A B) and the Coulomb energy; the grid sizes are 
determined by the hole density alone. These results appear to be consistent with the DMRG numerical solution of 
the t — J modeled which also finds loops of holes, except that in our case the periodic boundary conditions and the 
Coulomb interaction yield a "tile grid" as opposed to "droplets." Note the almost perfect (infinite charge correlation 
length, ^c) crystal structure obtained (Fig. @c). It is noteworthy that a typical solution yields a finite dipole moment at 
each grid intersection, which, in turn, can take one of the two orientations (along two diagonals of the computational 
box), thus creating a highly degenerate system of moments (see Fig. ^ below). 

We recall that the presented solution is obtained assuming an arbitrary dipole orientation with a constant hole 
dipole magnitude, i.e., a continuum of angles between the dipoles and a fixed axis. As explained in the previous 



section, at very low doping the dipoles would reside near the BZ diagonals, i.e., they would assume almost "discrete" 
orientations. In order to study the effect of this "discreteness" we have performed the same simulation this time 
assuming that hole dipoles can take only one of four directions (determined by the momenta of the maxima of the 
lower Hubbard band). We find that the physics of the pattern formation is qualitativelly unaltered (hence we do not 
present them here), except for one important difference: the "bending" of stripes at the grid intersections disappears, 
i.e., one no longer has a finite dipole moment at these intersections. 

Another way of quantifying this ordered phase is by straightforwardly calculating the "string tension," which, at 
T ^ is equal to d^U /dx^ , where U is the total potential energy and a; is a small hole (or stripe) displacement; a 
large string tension indicates a high stability of the obtained phase, and vice versa. In Fig. ^ we show the string 
tension for motion perpendicular to a grid side compared with the motion along a side. As seen in the figure, the grid 
phase (and, as discussed below, the stripe phase) is extremely stable with respect to the hole motion perpendicular 
to the holes line segments, due to the Coulomb interacton. On the other hand, at larger doping values and fixed hole 
density the stripes are almost compressible, i.e., the motion of holes along a stripe is rather soft. The anisotropy of 
the perpendicular and longitudinal string tension decreases with decreasing B. 
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FIG. 4. Top: Energy as a function of a hole position, reflecting the string tension in the stripe, for B = A eV and n = 15%. 
Clearly, the motion of holes perpendicular to the the grid directions is quite hard and while the motion along a line stripe 
is much softer (see also text). Bottom: average string tension as a function of _B, at n = 15% Note the almost exponential 
dependence (solid line) up to the critical value of B ~ 4 eV where the system undergoes the first order transition between the 
ferro- dipolar and grid phases. 



At fixed AF correlation length the four observed phases yield a diagram which we show in Fig. We remark 
that in all phases a non-vanishing value of A leads to a decrease in the effective value of B at which the transitions 
occur, as shown in Fig. The isotropic term A alone newer produces any non-trivial geometric phase (e.g., stripes), 
even with inclusion of lattice effects. We find that the transition between the ferro-dipolar and the grid phases is first 
order, as indicated by the coexistence of phases in Fig. |^. Note that this transition always occurs on increasing the 
doping level to sufficiently (and artificially) high values, where our theory need not apply. No coexistence of phases 
has been observed at other transitions, suggesting that they are second order. We also recall that our calculations 
are quasiclassical and thus the obtained geometric (stripe) phases are insulating. Moreover, in our formalism the hole 
density within a stripe can assume an arbitrary value, depending on the dipolar interaction strength. 
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Short Range Attraction A (B=2 eV) 
FIG. 5. (sl) A — 0, n — B phase diagram (from Ref. onlinecitesybcj). (b) Fixed B, n — A phase diagram. It is noteworthy 
that the experimentally relevant values of A and B are of order ~ 1 eV. The slightly higher values of B, at which we find 
geometric (grid and stripe) phases, are due to the fact that we consider unscreened Coulomb interaction and in reality they 
would be considerably smaller (see also Fig. 




FIG. 6. Coexistence of dipolar and grid phases, indicating the first order nature of the transition between them. 

In the cases presented above we have assumed a uniform magnetic dipolar interaction. It is well known that there are 
orthorhombic and tetragonal distortions in practically all transition metal oxides. In particular static stripe formation 
has only been observed in the low temperature tetragonal phase of Lai.6-2:Ndo.4Sra;Cu04B. In order to study the 
influence of the anisotropy we assume that the magnetic dipole sizes along the x and y directions have anisotropy a 
(a = 1 corresponds to the isotropic case). Figure ^ shows the pattern we obtain for a = 0.8: 
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FIG. 7. With a small x — y directional anisotropy in the dipolar interaction, Eq. (y), with the anisotropy parameter a = 0.8 
(in the uniform case a = 1), the holes (shown by circles) form a stripe, rather than a grid pattern. Note the hole dipole 
orientations (shown by the line segments), altering direction in neighbouring stripes, corresponding to the n phase shift of the 
local magnetization between magnetic domains separated by stripesE3 



The rotational symmetry is broken and a stripe superlattice is formed, with a charge ordering vector K = (27r/£)x, 
where £ is the inter-stripe distance. More importantly, the total dipole moment in this state vanishes. As explained 
in Sec. ||, this yields a Fourier transform of the magnetization S{q) — (5*2 (q)) peaked at Q ± K in momentum space. 
Of course, it is reasonable to assume that in twinned single crystals, used in inelastic neutron scattering experiments, 
one has domains which average out this anisotropy. Note that in both calculated geometric phases (see Figs, js] and 
0), the interstripe distance is much larger than the intrastripe distance, in agreement with experimental findings in 
underdoped cuprates. 

Our results are somewhat sensitive to the applied boundary conditions. First, for a small computational box the 
exact size of the grid depends on its commensuration with the box length, which, in turn depends on the density. On 
increasing of the size of the computational box, the grid size depends only on the physical parameters, as explained 
below Fig. ^. In addition, for a large computaitonal box the grid pattern, shown in Fig. |^a, acquires point or line 
defects, shown in Fig. 
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0.25 



FIG. 8. (a) A typical low energy state in a larger computational box, obtained for the same value of B and n as in Fig. ^. 
Note the presence of point and line defects. It is important to realize that the average size of the grid units depends only on 
the physical parameters and very weakly on the size of the computational box. (b) The charge density, pk, averaged over a 
number (of order 10) of hole configurations. Note that the peak positions are the same as those shown in Fig. ^ but due to 
the presence of defects the intensity of the higher Bragg peaks vanishes. 

This leads to the reduction in the higher order peaks observed in Fig. with no change in their wave-numbers, 
indicating the finite value of as seen in Fig. |^. The sensitivity to boundary conditions is further seen in finite size 
calculations, i.e, not assuming periodic boundary conditions, but with an appropriate neutralizing charge background. 
In this case the holes do not form geometric phases, although they still form stripes, as seen in Fig. |a. However, in 
a finite system even very small anisotropy (a ~ 0.95) again leads to stripe formation, as seen in Fig. |b. It is worth 
mentioning that the stripe formation occurs for much smaller values of B and the same density and AF correlation 
length in a finite system. 
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FIG. 9. (a) A typical low T metastable state of the hole system, obtained without periodic boundary conditions and with 
an appropriate positive background included. While the ground state is still a geometrically ordered state, it is practically 
unreachable due to the presence of many metastable states, (b) Adding very small directional anisotropy (a = 0.9, with a = 1 
corresponding to the uniform case shown in panel (a); see also Fig. ^) yields stripes with some defects. Panels (c) and (d) show 
Pfe as a function of momentum corresponding to the results shown in panels (a) and (b) , respectively. 

We have also performed simulations in the presence of a realistic underlying periodic lattice and have found that 
this creates slight distortions in the phases, pinning loops more strongly. In particular, the peaks in /9(q) sharpen in 
some of these phases. 



A. Role of disorder 



We now turn to the effects of point disorder. There are several kinds of impurities which are of experimental interest 
in transition metal oxides. We divide them into four distinct groups: 
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FIG. 10. The effects of impurities of different types. The four panels show hole (circles) and impurity (stars) positions, 
for the case of: charged impurities, placed d = 6Aout-of-plane (top left), in-plane repulsive uncharged impurities (top right), 
in-plane repulsive charged impurities (bottom left), and in-plane impurities with a local magnetic moment which destroys local 
AF ordering (bottom right). Clearly, all impurities destroy the stripe order, although the out-of-plane impurities and the 
uncharged in-plane impurities are nearly as effective as the in-plane charged impurities. The latter lead to the a glassy phase 
of stripe segments at relatively low concentrations, of order a few percent. 



(I)_Qut-of-plane impurities, such as Sr in LSCO compounds, (II) in-plane charged impurities, such as (presumably) 
LijLj (III) in-plane uncharged impurities, such as Ni and (IV) in plane uncharged impurities which induce a magnetic 
moment, such as Zn. Hence, we model the impurity effects by adding a random (in position) potential to the 
Hamiltonian (Q), which is either short-ranged and located in plane or, as in the case of charged impurities, long- 
ranged (Coulomb) and either in-plane or out-of-plane (a distance d from the plane where the holes are located). In 
the case of type IV we have also altered the dipolar interaction in the vicinity of an impurity, i.e., the magnetic 
interaction is multiplied by a factor tanh(r/i?i), where r is the distance of a hole to a nearby impurity and Ri is the 
effective radius of the. impurity, which, for the case of Zn, has been estimated to be of order 2 lattice sites around 
each impurity atom.E^ Examples of the effects of the four types of impurities are presented in Fig. |l^. 

Clearly, all four types of impurities lead to the destruction of the geometric (stripe or grid) hole order at sufhcicntly 
large impurity concentration, Ci . On the other hand, in all four cases stripe ordering persists through the formation 
of line segments of holes, resulting in a new, glassy phase.E3 Moreover, the four impurity types exhibit different mech- 
anisms for destroying the stripe order. The charge ordered phases are practically unaffected by a small concentration 
of uncharged impurities (see Fig. p^), i.e., the stripes simply avoid impurity sites. Consequently, the stripes persist 
to relatively high concentrations of this type of disorder. 

The charged impurities first lead to stripe deformation, i.e., the stripes pass either very close to the impurities (for 
attractive, pinning impurities) or very far from the impurities (for repulsive impurities), in order to maximize the 
potential energy (see Figs. [lO^ and p^). With increasing Ci the stripes rupture and only stripe segments persist. 
Finally, the impurities with a local magnetic moment affect the formation of the spiral spin phase, responsible for 
the (attractive) dipolar interaction. Since the magnetic interaction is strongly suppressed in the vicinity of such an 
impurity site, even the stripe segments cannot exist there, as shown in Fig. [lO|d. 

Impurities are especially effective in destroying the ordered phases found at small B. For example, the Wigner 
crystal state becomes glassy at relatively low impurity concentrations. This happens because, e.g., in the case of 
impurity type I, the attractive Coulomb energy between impurities and holes scales like e^/d, where d is the distance 
between the planes in which the impurities and holes reside, while the average inter-hole Coulomb energy behaves 
like e^^fal- Thus when Ug < 1/d^ the holes are pinned by impurities. 

In general the role of impurities depends strongly on the impurity concentration, a. However, the magnetic dipole 
interaction is sufficient to retain the main orientation, as seen in Fig. ^ where we have plotted the correlation length 
as a function of c^. This leads us to conjecture that with the addition of the kinetic energy the holes can move in 
string segments in an orientation given basically by the phase diagram of the clean system. The stripe motion would 
then be caused by mesoscopic thermal or quantum tunneling of the finite strings between the jainima of the overall 
potential. This would lead to non-linear field dependence in the low temperature conductivity.EJ 
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FIG. 11. Zero temperature charge correlation length as a function of impurity concentration for the in-plane charged impurity 
case: the result is obtained by measuring the static correlation length, and then averaging over many (of order 20) different 
impurity configurations. 



B. Finite temerature results 



We now proceed to the finite T results. The numerical procedure is identical, except that the temperature is lowered 
adiabatically to a finite value (i.e., a numerical annealing). In a classical simulation this is equivalent to introducing 
kinetic energy into the system. 

In the case of a Wigner crystal, i? = 0, we find that the introduction of finite T melts the crystalline structure and 
the resulting phase is the 2D Coulomb gas. The diagonal (glassy) phase is also unstable at relatively low temperatures. 
On the other hand, the geometricaly ordered states, for as <C where as is the distance between holes within a stripe 
segment, are all stable up to T of order ~ B/asoj.. At even higher temperatures, the stripe array melts with a 
temporal intermittency of the observed pattern: i.e., spatio-temporal intermittency. Fig. y_2| shows four stages of 
this melting process. We observed that the stripe melts through a rupture which results in creation of finite stripe 
segments that eventually (at constant and high T) disperse into individual holes. 

Note that the temporal geometric pattern (panel (b)) is not the same as that in the ground state. As mentioned 
before, there are many low lying geometric states, close in energy to the ground state, which can temporarily occur 
at finite T. Hence the dynamics of the stripe ordering is similar to that observed in glasses, characterized by non- 
gaussian fiuctuations. To show this we follow the dynamics of the hole system at temperatures slightly helow the 
melting temperature: we start from a low lying metastable state, such as that depicted in Fig. |l2|b, increase the 
temperature adiabatically to the point at which the structure begins to melt (which is a measure of the activation 
energy) and let the system equillibrate. 
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FIG. 12. Spatio-temporal intermittent behavior of liole ordering near the melting transition: three snapshots are shown. 
Panel (a) corresponds to a state which is clearly withing the basin of attraction of the ground state (the same pattern, albeit 
deformed), while panel (b) shows a state which is within the basin of attraction of another low lying geometric state with 
more dense stripes along one of the axis. Panel (c) shows the melted nematic crystal-like phase with the hole dipole moments 
aligned. 



In Fig. |l3|a we plot an energy histogram at this temperature (with the energy shifted by an arbitrary additive constant) , 
thus indicating an intensity of the energy states (bands): 
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FIG. 13. Top: a histogram of the energies at constant (high) temperature of order of the activation energy for a transition 
between the two energy basins, depicted by the states shown in Fig. The left maximum corresponds to the ground state, 
Fig. [l^, and the right to a family of low lying excited states. Fig. [l2p). Neither peak can be fit by a simple gaussian, indicating 
a glassy nature of the ordered states. Bottom: The average "power spectrum" corresponding to the result shown in the top 
panel: the vertical axis shows the sum of the squares of the Fourier components of the potential energy, within the nth octave 
(components 2"~^ through 2" — 1, for n > 0), averaged owr many (of order 300) spectra. The flat low frequency behavior 
(up to about the 9th component) is very close to a genericEj 1// noise and corresponds to slow fluctuations involving many 
particles, such as those yielding the transitions between the states depicted in Figs. [L2p, and |l2p. 



Obviously, there is a band of energy states, not far (fraction of an eV per particle) from the ground state, which 
are close in energy and metastable. These states are separated by a high barrier (the maxiniumjaf which would fall 
beyond the right edge on the plot), yet are close in energy, suggesting a rugged energy landscape.Lj Indeed, as shown 
earlier, formation of a string of holes creates a barrier for adding more holes to the string (they can be only added 
to the string ends). Thus, any geometrically ordered state (say, those with denser intra-stripe hole concentration and 
larger interstripe distances) must be separated by an energy barrier from other geometrically ordered states and in 
particular from the ground state. 

The potential energy states obtained suggest that the dynamics of stripe motion should be strongly governed by 
these low lying states and thus show a non-trivial fluctuation spectrum. Indeed, in Fig. ^3|b we show the power 
spectrum of the energy fluctuations for the solution described by the histogram in Fig. and see that the noise 
spectrum contains a strong 1// component for approximately two and half decades of frequencies. This indicates slow 
fluctuations, which we ascribe to collective motions of melted hole string segments. 

Another way of characterizing the melting of stripes is by counting "free holes:" in Fig. ^wc show the percentage 
of holes which are not in a part of an ordered pattern, as a function of T. As one can see, at the transition point only 
a small fraction of holes does not belong to a string segment, in agreement with our observation that the stripes melt 
by rupturing into smaller segments. 
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FIG. 14. Percentage of holes which are not a part of a stripe segment, as a function of temperature, for B — 4eV. 

Further study of the glassy dynamics of charge ordered phases in terms of the (free) energy landscape will be 
presented elsewhere. 
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IV. CONCLUSIONS 



In summary, on employing the SDW picture of transition metal oxides, we have studied the short-range and dipolar 
attractive forces generated by the AF fluctuations, together with long-range Coulomb forces. We have developed a 
novel numerical technique, which enables us to treat doped hole vacancies at finite concentration. We have studied the 
competition between long-range and short-range interactions and its influence on hole ordering in layered transition 
metal oxides. We have found a rich phase diagram for the clean system which includes a Wigner solid, diagonal 
stripes, grid (loops) and a macroscopic phase separation. For intermediate values of magnetic interaction this phase 
diagram is consistent with several different experimental measurements, such as the inelastic neutron scattering. In 
addition, on adding a small, but finite amount of anisotropy to the dipolar interaction we find that the ground state of 
the system of holes is the striped phase, found in La2-y-a;NdySr2,Cu04. In the geometric phases with strong magnetic 
interaction strength we have found a large string tension for the motion of holes perpendicular to the stripe direction. 
This is due to the Coulomb interaction and indicates strong stability of the obtained phases. 

We have also found the system of holes to be quite sensitive to the presence of charged impurities. In particular, 
adding out-of-plane attractive impurities pins the holes and, for small pinning energies, increases the melting tem- 
perature of the stripe phase, although it does yield a finite charge correlation length. In general, charged impurities 
are very effective in destroying the stripe order, especially those residing in the same plane as the holes, regardless of 
whether they are attractive or repulsive, although the stripe phases survive as finite stripe segments up to relatively 
high impurity concentrations. This suggests that nonlinear conductivity should be prevalent. 

The resulting hole patterns arc the result of frustration (competition between short-range and long-range forces): 
this frustration leads to collective motions, involving large number of particles which ultimately lead to geometricly 
ordered ground states. We have also studies the dynamics of the geometric phase formation and its melting. We find 
that the dynamics is characterized by a "glassy" behavior in that the energy landscape is rugged, as characterized by 
the spatio-temporal intermittency of the observed behavior. More importantly, we find that, for fixed (large) size of 
the magnetic (dipolar) interaction, there are fewer number of sharper minima, while the string tension of the stripes 
is larger. As a consequence, in this case the melting of the stripe phase occurs at higher temperature with increased 
doping concentration. 

The energy landscape is also charactecized by formation of domains, separated by defects. This picture is in 
agreement with recent NMR experimentsllil in which small activation energies are easily attributed to domain growth 
and/or motion. Thus, further study of our model will include the dynamics of the domain growth and their melting. 
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APPENDIX A: MAGNETIC DIPOLES OF THE HUBBARD MODEL 

In this appendix we study the Hubbard model, Eq. (|l|) in the SDW state. Our approach is similar to that presented 
in Refs. 28 3^. Hence we only briefly review the calculation leading to Eq. (|5|), the central equation of the paper. 
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The relevant order parameter in our case is the spin density in the z direction: 

5^(q)=5]a4+q_„Ck.o, (Al) 

k,a 

which in the SDW state has a finite expectation value at q = Q = (7r/a, 7r/a) because of the nesting of the half-filled 
Fermi surface. In this case the mean field Hamiltonian reads 

HmF = e/c4_„Ck,a - Yl "4+Q,aCk,a , (A2) 

k.a k,CK 

where 

5 = 1(S^(Q)), 

Efc — —2t {cos{kxa) + cos(fcya)) . (A3) 
The Hamiltonian ( |A2| ) can be diagonalized immediately using the Bogoliubov transformation: 

7k,a = «kCk,a - aWkCk+Q,Q (A4) 

where 



Vk = 



A = -— . (A5) 



In this case the mean-field Hamiltonian reads 

Hmf =Y.^^ (7k!a7£,o - 7k!o7^,a) , (A6) 

k.Q 

where the sum over k is restricted to the magnetic Brillouin zone and, at half filling, the ground state |0) is defined 
such that 

7£.J0>=0 

7k!jO)-0. (A7) 

Thus, at half-filling the conduction band is empty and the valence band is separated from it by energy A which is 
the Mott-Hubbard gap. It is known that this theory recovers the results of the Heisenberg model very well. Consider, 
for instance, the average spin density in ( |A1[ ) in terms of the new operators (recall that the conduction band is empty 
and therefore does not contribute) 



(5^(q)) = -2^Uk+q-Q^k(7klq-Q,.7i^,o> 



k.Q: 



k k 



which, together with (A5), yields the gap equation: 



1 1 1 
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Since we are going to consider hole doping we can neglect the terms involving the conduction band operators for 
temperatures T <C A. As shown in ^ new interactions are generated by the an tiferromagnet in the presence of the 
holes, given by and Hxy The non-interaction hole Hamiltonian is, from (A6): 

Ho = -Y.E>'<c.llc.- (AlO) 

k,Q 

Close to the half-filled Fermi surface one sees that the hole mass is 

For any state j^P) of the system we can define the hole operators as 
in which case the Hamiltonian reads 

iJo = 5^£;fc<„Va (A13) 

plus unimportant constants. 

The interacting parts of the Hamiltonian can also be written in terms of this new operators. For instance. 



= [(^-(Q) - K(2k- Q)/4)m2 , - V,{2k)llJA] hl ji^,^ 

k,Q: 

k,k' 

+ V,{k-k' + Q)m2 „/i^/ik',a/iLk.^/i-k',/3) , (A14) 



where the sum over spin indices is implicit and 

VM = JI^2(M_ (A15) 

with 

Xo(9,'^) = -^E(i- EkEk+g ) {u - E,,+g - Ek ' u + Ek+g + Ek) ' ^^^^^ 

A similar expression is valid for the transverse components of the interaction 

4 ■ 

N ■ 



- E - ^>lk' - V+- {k-k' + Q)pI,,] 



k,k' 



X 4,a'^Ja''»k%a'/lLk,/3Cr/3,/3'/l-k%/3' , (A17) 



where is given by an expression similar to (A15) with Xo replaced by 



==-^y(l- ''''^^~^" )( 1 ] . (A18) 



2N^\ EkEk+q J\uj-Ek+q-Ek UJ + Ek+q+Ek 

Moreover, the coefficients that appear in these expressions are defined by 

nT-fc,fc' = UkVk> + VkUk' 

lk,k' = UkUk' + VkVk' 
Pk,k' = UkVk' - VkUk' 

nk,k' = UkUk' - VkVk' ■ (A19) 



21 



Observe that the interactions renormahze the dispersion of the holes as well, that is, Ek E^. While is 
essentially a short range attractive interaction in the spin-symmetric channel, H^y has two components: one of them 
is also an attractive interaction in the spin-symmetric channel but the other is a long range dipolar interaction which 
depends on the momentum of the hole. 

The Hamiltonian, which consists of terms given by Eqs. (A14) and (A17), is very hard to deal with. One notices, 
however, that the mean field energy Ek is degenerate along the magnetic Brillouin zone. This is an artifact of the 
theory and the degeneracy is broken by any small perturbation such as the corrections discussed before or a next 
nearest hopping, i', for instance. In this case, the dispersion has a minimum at (±7r/2, ±7r/2). Thus, in order to 
study the long wavelength limit of 
paper by Schrieffer, Wen and Zhan; 



e theory it is sufficient to focus on these points of the Brillouin zone. In the 
the authors focused entirely on the part ot-the Hamiltonian since the form 
factors rt/c.fc' and Pk,k' vanish at the Brillouin zone. As shown by Frenkel and IIankc,c2l if one keeps the leading order 
in momentum we can write 



Pk,k' 

and therefore the interaction term becomes 



^+-(q+Q)pL+««2C/ 



1 1 



K) + (ky-K)\ 



= 2U 1 



.qxQy 



(A20) 



(A21) 



which has dipolar form. 

Thus, in the first quantized language the interactions have the form 



Hr 



[Aa-(ri)(7^(r2) - B (a+(ri)a- (ra) + a-(ri)a+(r2))] ,5(ri - ra) 



(A22) 



where a are spin operators. Notice that the singlet state | t: i) ^ I i: T)j clearly minimizes the energy of interaction 
between the two holes. In this case we end up with the charge interactions only. 

From our simulations we see that the charge orders with some characteristic vector K such that 



(A23) 



acquires a finite expectation value at q = K. Thus, we can always write down a mean field version of (A22) plus the 
long range Coulomb interaction as 



Hi 



4 



^'k.a 1 



(A24) 



where Vk has to be calculated from ( [A14| ) and ( |A17| ) and 



Observe that in this case the Brillouin zone is further reduced and we can define new operators 



with 



Wk 



tk 



+ tkhl 

tkhk, 



Wkh) 



k+K,Q 

t 

k+K,Q 



(A25) 



(A26) 



Ei = 




(A27) 
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and the Hamiltonian is diagonal 



H = 

k,a 



E - , (A28) 



where the sum is done in the new Brillouin zone. Self-consistency requires that 



= ^ E ^'^^^ («l<a> - (^k Vk, j) (A29) 



P N 



which can be evaluated for any hole-filling. 

Now let us go back to the issue of magnetization which is important for neutron scattering. From (|A8|) one has 



(5^(q)) = -2E"k+q-Qt'k(7klq-Q,.7ii,a) 

= -2(5q,Q+K ^k+K^k(/^k+K,a/^k.a) 
k,Q: 

- 25q,Q+K E "k+KWk^fcWfc ({dt,o,d^L) - (A30) 



k,a 



and one sees that the magnetization is now peaked around Q + K instead of Q. For stripes aligned along the y 
direction this is possible of course when 



27r 

K = ±^x, (A31) 



where £ is the inter-stripe distance. 



APPENDIX B: DIPOLES OF THE T J MODEL 

Let us consider the SDW theory of the t — J model d la Shraiman and SiggiajE^I described by 

-ff - -t E + h.c. + J ^ S, • S, (Bl) 

(iJ) (id) 

and use the slave fermion representation 

Ci,a = '4^a,r^a.M,y , (B2) 

where V'l i creates a hole (fermion) on a site i in the sublattice a = A,B (which labels the " spin" of the hole) and 
Za,i,a is a Schwingcr boson on the same sublattice (spin wave). In order to obtaiu-the dynamics of the holes alone we 
trace out the spin-wave degrees of freedom. The static part of the interaction isEJ 

^ss^-]^ E nk,k',q)Vi(k)V'B(k + q)^l3(k'+q)VM(k'), (B3) 

k,k',q 

where the momentum sum is restricted to the magnetic Brillouin zone and 

(Ak — Ak+q) (Ak' — Ak'+q) 



nk,k',q) = -9 



1 - Aq 



(Ak + Ak+q) (Ak' + Ak'+q) 
^ 1 + Aq 

with 



(B4) 
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cos{qy)) . 



(B5) 



The coupling constant g is a function oLt and J. In the strong couphng Umit (t << J) 5 ~ 8<^/J, while in the weak 
coupling limit {t >> J) we have g w JEB. 

Observe that (B3) does not have a kinetic term for the holes. The kinetic energy has to be obtained from the hole 
self-energy at zero frequency and can be written as 



(B6) 



where Ck has a minimum at (±7r/2, ±7r/2j0. For a low density of holes these are the only points of interest and 
therefore we can look at the interaction (B4) strength close to these points. Observe that at these points we have 
Aq ^ and therefore the interactions are dominated by the first term in (B4) which describes the fluctuations of the 



staggered magnetization (with characteristic wave- vector Q — (tTjTt)). The second term describes the fluctuations 
of the homogeneous magnetization (q = 0) which is not of direct interest here. In this case, for k and k' close 
(±7r/2, ±7r/2) to the interaction can be approximated by 



y(k,k',q)«-g 



Ak+qAk'+q 
1 - An 



(B7) 



The problem can be further simplified if one works with the upper half-part of the original BZ instead of the magnetic 
BZ, as shown in FigJl^. This can be accomphshed by a shift of the lower part of the BZ by Q. 
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FIG. 15. Choice of the BZ. 
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Moreover, working in the long wave-length limit, that is, with q ^ 0, one sees that there are four values of V^(k, k', q) 
of relevance: 1) when k = k' = Q/2 



2) when k = k' = Q*/2 (where Q* = (-tt, tt)) and 



V'22(q) ~ -9 



3)when k = Q/2 and k' = Q*/2 



4)and finally k = Q*/2 and k' = Q/2 



^12 (q) « g 



^21 (q) -1^12 (q) 



(B8) 



(B9) 



(BIO) 



(Bll) 



We can now split the sums in (B3) to the regions around these two points and introduce a cut-off in the momentum 
sum A such that q << A << tt. In this case the Hilbert space of the problem is divided into two different sub-Hilbert 
spaces and the hole operator can be rewritten as 



COS ( ^ • r ) -I- ■0Q,i,2 COS 



Q* 



(B12) 



where 



q 

V'c.,,,2«$]e'''--'Vo(^ + q). 
q 

It is convenient to define the operator 

= II'/'ka(k + q)V'Aa(k) 

k 

Pl(q) - E AaO^)Ma{k + q) = E Aai^ ' 



(B13) 



(B14) 



with a = I, 2. Using these new operators and results (B8)-(Bf f ), we can rewrite the Hamiltonian ( p33| ) as 

^^^^ ^ E ^ [i' (pl(q)pi(q) +P2(q)p2(q)) 

q 

- ill - qD (pI(q)P2(q) +pj(q)pi(q)) + 2q,qy (pj (q)pi(q) - 4(q)p2(q))] ■ (B15) 

This does not have a very transparent form. In order to see that this Hamiltonian has the form of a dipole-dipole 
interaction we define the vector operator 



D (q) = {pi (q) - P2 (q) , pi (q) + P2 (q) ) 
Dt(q) = -L (pU-q) -4(-q),pl(-q) +4(-q; 



(BI6) 
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and rewrite ( B15| ) as 



q 



Dt(-q).D(q)-2 



(Dt(-q).q) (D(q).q) 



(B17) 



Observe that the first term in ( B17 ) is q-independent and will lead to a local interaction which has the usual scalar 
form. We are interested in the second part of the Hamiltonian. Defining the Fourier transform 



Z?(r) = ^e''i"-i?(q) 



(B18) 



the second term in (B17) acquires the required form 



Hdd =^ g dr dr' 



D+(r)-D(r')-2 



(Dt(r) -r) (D(r') -r')' 
(r — r')2 



(B19) 



which is the second quantized form of the magnetic dipole-dipole interaction. 

In order to ppt this Hamiltonian in a form involving the magnetic dipole defined by Shraimann and Siggia consider 
their definitionE3: 



(B20) 



k,a,6 



where fj, — x,y, z refers to indices of the Pauli matrices t^j, and therefore act in the sub-space of the sublatticcs A 
and B and a — x,y refers to the space indices. In particular, the lowering and raising operators associated with this 
magnetic dipole operators have the interesting form 



P_,c«(q) = 2^ sin(fc„)V'l5(k + q)V'A(k) 

k 

(q) = 2^sin(fcQ)'0]i(k + q)'0s(k). 



(B21) 



In the approximation we are employing we can split the summation in ( |B21 ) around Q and Q* in order to get (this 
can be done because the sine function is smooth around these two points) 





.2J2[AbA^' 

k 


hq)VM4(k) 


~AbA^^ 


hq)VM,2(k) 


P-A<i) ^ 


.2^[vL,i(k^ 

k 


hq)VM4(k) 


+ AbA^^ 


hq)VM,2(k) 


p+A^) - 


k 


-q)V'B,i(k) 




-q)^/'B,2(k) 


P+A(i) ^ 




-q)^'B.i(k) 


+ AaA^^ 


-q)V'i3,2(k) 



(B22) 



where the signs come from value of the sines around the two points in the FS. These dipole operators can be also 
written trivially in terms of the operators in (B14): 



P_,,(q) =2(pi(q)-p2(q)) 
P_,,(q) -2(pi(q)+p2(q)) 



(B23) 



and so on. By direct comparison with (|B16[ ) one finds 



D(q) = ^(P-,.(q),P-,,(q)) , 



(B24) 
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which makes clear the connection. On taking the Fourier transform of the magnetic dipole operators back to real 
space one finds, for instance, 

<,2,ii^A,2,^) ■ (B25) 

This cxplicitcly justifies our earlier claim that the dipolar interaction is due to the coherent hoping of holes between 
two different sublattices (at the same position in space). 
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